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We describe a novel method to obtain thermodynamic properties of quantum systems using 
Baysian Inference - Maximum Entropy techniques. The method is applicable to energy values 
sampled at a discrete set of temperatures from Quantum Monte Carlo Simulations. The internal 
energy and the specific heat of the system are easily obtained as are errorbars on these quantities. 
The entropy and the free energy are also obtainable. No assumptions as to the specific functional 
form of the energy are made. The use of a priori information, such as a sum rule on the entropy, is 
built into the method. As a non-trivial example of the method, we obtain the specific heat of the 
three-dimensional Periodic Anderson Model. 
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The problem of obtaining thermodynamic properties 
from Quantum Monte Carlo (QMC) simulations is one 
of long-standing interest. ||]. Although the internal en- 
ergy, i.e., the expectation value of the Hamiltonian, is 
one of the easiest quantities to obtain via QMC, the free 
energy is almost impossible to obtain directly in a simu- 
lation. Likewise, the specific heat, i.e., the temperature 
derivative of the internal energy, is very difficult to obtain 
directly. Hence, one must turn to indirect methods. 

Several methods to obtain the thermodynamic proper- 
ties of model systems via QMC have been proposed, but 
all suffer from limitations of one sort or another. To a 
large extent, these stem from the use of a specific func- 
tional form to fit the internal energy of the system. We 
propose a novel method to obtain the internal energy, the 
specific heat, the entropy, and the free energy as a func- 
tion of temperature via QMC which does not impose any 
functional form on these quantities and alleviates several 
other problems in the current methods. Our technique 
relies on probability theory and Maximum Entropy to 
obtain the most probable thermodynamic functions con- 
sistent with the QMC data and prior knowledge, such as 
a sum rule on the system's entropy. 

In the remainder of this paper, Sec. O reviews cur- 
rently used techniques to obtain thermodynamic prop- 
erties from QMC, their limitations, and the desirable 
features of a new technique. Sec. Ill contains a brief 



overview of our method to obtain thermodynamic quan- 



tities from QMC data. In Sec. IV, we review the theo- 



retical underpinnings of the method. Maximum Entropy. 
Sec. M sets forth the algorithmic details of our method. 
In order to test our method, we apply it to a non-trivial 
problem - the 3d Periodic Anderson Model ~ in Sec. yf\. 



Our summary is given in Sec. VII 



The free energy and its derivatives, including the spe- 
cific heat, provide experimentally relevant insight into 
a system's temperature evolution and phase transitions. 
Unfortunately, both direct and indirect QMC measure- 
ments of such quantities are notoriously difficult to make. 
To understand why, we discuss two methods for obtaining 
thermodynamic quantities in this section. 

Typically, one tries to obtain the thermodynamic prop- 
erties of a system by performing QMC simulations at 
various discrete temperatures, then fitting the resultant 
energy data to a functional form. Generally, this func- 
tional form is not known, so a physically-motivated form 
must be chosen. The recipe is to fit the internal energy 
E{T) to a functional form, which may then be differen- 
tiated explicitly to obtain the specific heat 



C{T) 



dEjT) 
dT ' 



(1) 



From the specific heat, the entropy S(T) may be calcu- 
lated by integrating 



S{T)= f dT 
Jo 



, C{T') 
T ' 



(2) 



Then, the free energy F may be obtained from the rela- 
tion 



F{T) ^E{T)-TS{T). 



(3) 



While apparently sound in principle, this prescription 
can manifest several serious problems. For one, the 
derivative in Eq. 1 enhances the statistical uncertainty 
in the fit. At low temperatures the procedure is further 
complicated by the division by T in Eq. 2. Similar prob- 
lems emerge at low T when the specific heat is evaluated 
directly C = Ue'^) - {Ef) /T^. Thus, there is no guar- 
antee that the total entropy obtained by this method 
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will equal the total infinite-temperature entropy of the 
system. If not, then both the specific heat C{T) and the 
free-energy F{T) may be unreliable. 

One technique using this prescription is to fit QMC in- 
ternal energy data to a pair of functional forms, one for 
low-temperature data and another for high-temperature 
data. Q This method has been successfully apphed, but 
it does have several important drawbacks. First, separate 
functional forms in the nature of polynomials are used 
for the low- and the high-temperature regions. These 
are chosen after viewing the data, based on the shape 
of the E{T) curve. Presumably, care must be taken to 
avoid spurious features in the derivative C(T) of the in- 
ternal energy E(T) at the point where the two polyno- 
mials are joined. While an analytic function such as a 
polynomial expansion is a reasonable choice for an inter- 
nal energy function on a finite-dimensional lattice, there 
is no guarantee that sufficient terms have been chosen for 
the polynomial and the most reliable test is a goodness 
of fit. Consequently, this technique requires a great deal 
of costly QMC data, especially in the neighborhood of 
a phase transition, in order to insure that a reasonable 
goodness of fit is obtained. 

Another recently-proposed method to obtain thermo- 
dynamic properties from QMC data is to fit the internal 
energy E{T) to a physically-motivated functional form. 
P,H In this method, appropriate for a lattice simulation, 
the internal energy is fit to a sum of exponentials 



EiT) = Eo + Y, Cut 



■lA/T 



(5) 



The specific heat and the entropy are then obtained ac- 
cording to Eqs. 2 & 0, respectively. 

One way to view the exponential functional form of 
this method is to note that physically, one may expect 
different energy scales to become important as the tem- 
perature T is varied. At those energy scales, contribu- 
tions to the internal energy E{T) are effectively switched 
on. This is manifest in the fitting parameter A; at each 
temperature corresponding to the energy scale nA for 
each of the n terms in the expansion, another term in 
the expansion contributes to the energy. The amplitude 
of this contribution is set by the related coefficient c„ . 

This technique has been successfully applied. It has at 
least one major advantages over the polynomial method 
in that it does not splice together two functions for differ- 
ent temperature regimes. Nevertheless, it also relics on 
a goodness of fit test to determine whether a reasonable 
number of fitting parameters have been chosen. Since it 
uses a gapped form for the internal energy, it is not suit- 
able for systems in the thermodynamic limit, as may be 
studied using dynamical mean field techniques. Further- 
more, both the polynomial and the exponential fitting 



schemes are inaccurate when the number of fitting pa- 
rameters is small and become ill posed as the number of 
fitting parameters become large. Thus it is difficult to 
determine how many coefficients to use. 

We now present a method that overcomes these draw- 
backs and possible pitfalls. Our method incorporates ad- 
ditional a priori information that the energy E{T) in- 
creases monotonically with the temperature T so that 
the specific heat is positive definite, C{T) > 0, and in 
the form of the infinite temperature entropy 5*00- Thus 
it requires less QMC data than either of the functional 
fitting methods. It performs a search for the most prob- 
able energy E{T) given the data and prior information, 
and therefore removes the question of determining the 
number of fitting parameters appearing when using the 
other methods. It is applicable to both finite-dimensional 
lattice QMC as well as QMC from non-gapped systems. 



III. OVERVIEW OF THE NEW TECHNIQUE 

We start with the observation that given the appropri- 
ate distribution function K{P, co) relating the energy tv to 
the inverse temperature (3 — 1/T, one can write the inter- 
nal energy as a weighted integral with a positive-definite 
weight p{i-u) 



E{T) = - <kjuoK{f3,uj)p{uj). 



The specific heat is obtained by differentiating 



C{T) = 



dE{T) 
dT 



dK{P,L0) 

— df — ^^^'' 



(6) 



(7) 



One may then use Eqs. y and|3|to obtain the entropy and 
the free energy, respectively. 

The entire problem of obtaining the thermodynamic 
properties of the system then reduces to that of numeri- 
cally inverting the integral equation Eq. |6| for the weight 
p{Ld) given noisy QMC data for the internal energy E{T). 
This is a well-known problem for which there exists a 
well-developed, powerful technique, the Maximum En- 
tropy Method (MEM). |] (The MEM is discussed in the 
next section.) The kernel K((3,uj) corresponds to a blur- 
ring function which acts on a spectrum p(cj). Typical 
blurring functions for a quantum system are the Bose- 
Einstein and the Fermi distribution functions, as dis- 
cussed in detail in Sec. 



[V 



To understand the fundamental difference between 
other methods and our MEM technique, it is important 
to understand the questions which the two methods an- 
swer. The other methods rely on fitting noisy data to a 
functional form. They start with a physically-motivated 
functional form for E{T) and seek to find the most likely 
curve of the infinitely many curves which optimize some 
likelihood function such as x^. In reality, the functional 



form of the energy is not known. What is known is that 
finite temperatures spread the excitation spectrum of a 
quantum system. The question one might hke to ask, 
"what is the curve that fits the data" is therefore ill- 
defined. Without additional regularization, an infinite 
number of curves fit the data and, unless the precise func- 
tional form of the internal energy E{T) is known, there 
is no precise answer to this question. However, if we 
know the form of the thermal blurring function K{(3, uj), 
we may ask the question, "given the blurring function 
and any other relevant prior information that we know, 
what is the most probable spectrum p{cu) from which this 
energy data E{T) might arise?" As discussed in detail 
below, this is the precise question that the MEM answers. 
In addition to relying on fundamental properties of 
quantum systems instead of a functional form, further 
benefits also accrue from employing our MEM technique. 
For example, since the MEM gives errorbars on inte- 
grated quantities, the uncertainties in both E(T) and 
C{T) are known when obtained via our MEM technique. 
Other advantages of our technique will be discussed be- 
low. Before discussing algorithmic details, it is useful and 
instructive to briefly review the MEM. 



IV. THE MAXIMUM ENTROPY METHOD 

The Maximum Entropy Method (MEM) is discussed 
in detail elsewhere. [pP-M| Here, we wish to review only 
as much of the MEM as is necessary to understand the 
new technique. 

The MEM is frequently used to analytically continue 
QMC imaginary-time Green Function data to real fre- 
quencies. H However, it is a general technique that is 
not limited to analytic continuation or to QMC-related 
problems. In fact, the MEM has a relatively long history 
as an image reconstruction technique in photography and 
dynamic light scattering problems. pVjj^ 

In such problems, the observed image is the result 
of Gaussian blurring of light transmitted from a source 
through a medium, such as the atmosphere. Hence, the 
functional form of the image is not known and the ques- 
tion of whether the observed data flts a specific functional 
form is ill-defined - there are an infinite number of curves 
that fit the data! Instead, the best that one can do is to 
seek the most probable image given the data. This is 
exactly what the MEM sets out to accomplish. 

This is done using Bayesian statistics. If there are 
two events, a and b, then by Bayes' theorem, the joint 
probability of these two events is 



P(a,6) = P{a\b)P{b) = P{b\a)P{a), 



(8) 



where P{a\b) is the conditional probability of a given b. 
The probabilities are normalized so that 



P{a) = I dbP{a,b) and f daP{a) = 1. 



(9) 



In our problem, we search for the spectrum p which 
maximizes the conditional probability of p given the data 

E, 



P{p\E)=P{E\p)P{p)/P{E). 



(10) 



Typically, one calls P{E\p) the likelihood function and 
P{p) the prior probability of p (or the prior). Since we 
work with one set of QMC data at a time, P{E) is a 
constant during this procedure and may be ignored. The 
prior and the likelihood functions require more thought, 
and are discussed in detail in Ref. [pi, here we present 
the salient results of that discussion. 

If the spectrum is positive-definite, we may think of it 
as a un-normalized probability density: 



dujp(uj) < oo. 



(11) 



Then by Skilling, |p| the prior probability is proportional 
to exp(ckS') where S is the entropy defined relative to 
some positive-definite function m{uj) 

/oo 
duj[p{uj) — m[ijj) — p[ijj) ln(/9(w)/m(w))] , (12) 
-oo 

and 



P{p) oc P(p|m, a) oc exp(aS'), 



(13) 



where m{uj) is the default model since in the absence 
of data p = m. Selection of the default model for this 
method is discussed in Sec. |V|. The other unknown quan- 
tity a is determined during the MEM to maximize the 
probability of the image p given the data. 

The likelihood function follows from the central limit 
theorem. If each of the measurements Ei^T {Et = E{T)) 
of the energy at a specific temperature T is independent, 
then in the limit of a large number of measurements Nd 
to determine each Et the distribution of the Et becomes 
Gaussian. The probability of measuring a particular Et 
is 



P{Et) = -=^exp[-^((i?T) -£;t)V2], 
V ZnaT "t 

with an error estimate given by 

^^Y.i^ET)-E,^Tf 



(Jt 



NdiNa - 1) 



and 



(i^.) = ^E^> 



(14) 



(15) 



(16) 



for the Nd measurements of Ei^T at temperature T. 



Then the Hkehhood function P{E\p) of measuring the 
set of E for a given image p is 



P{E\p) 



= -xV2 



where 



X 



'^ {ET-E.^K{cu,T)p{co)y 



(17) 



(18) 



and we have discretized the integral Eq. y. 

We are now in a position to perform the MEM and 
find the most probable image p given the data Et- We 
wish to maximize the joint probability of the image or 
weight p given the data E; the default model m; and the 
Lagrange multiplier a 



P{p\E, m, a) ex P{E\p)P{p\m, a) 
_ exp(«5 - xV2) 
ZsZl 



(19) 



where Zs and Zl are normalization factors, indepen- 
dent of the image. For a fixed a and the given data 
E, the most probable image /3(a) is the one that maxi- 
mizes Q ~ aS — x^/2- This may be found, for example, 
using Newton's method. 

The details of implementing a MEM code and finding 
a are given elsewhere. We will not repeat that presen- 
tation here. A MEM code written according to Ref. g] 
is recommended for performing the technique we discuss 
herein. ||ll| Having discussed the general MEM formal- 
ism, we now turn to the specific algorithmic details of 
our new technique. 



V. ALGORITHMIC DETAILS OF THE NEW 
TECHNIQUE 

We desire to express the internal energy of the sys- 
tem as an integral over a density of energy levels times 
a relation between energy and temperature, according to 
Eq. g. To that end, we make the following Ansatz 



EiT) 



dujuj[F{l3,uj)pF{uj) + B{p,uj)pB{io)], (20) 



where F and B arc the Fermi- and Bose-distribution 
functions, respectively 



F{p,uj) 

B{(i,L0) 



1 



1 



1 



cpU] 



(21) 



for (3 = 1/T (we have set the Boltzmann constant equal 
to unity ks = 1). 

This Ansatz corresponds roughly to describing the 
energetics of the system as consisting of separate lin- 
ear contributions from Fermi- and Bose-excitations and 



imposes the constraint that the corresponding energy 
E(T) increases monotonically with temperature T so 
that C{T) > 0. In addition, since the degeneracy of the 
ground state and the total number of accessible states is 
generally know, the infinite temperature entropy is gen- 
erally known. 6*00 may be obtained from 



•^oo — 



dT 



duJUj[—pF{uj) 



dB 
df 



PBiu^)] (22) 



by noting that the temperature integral for the Fermi 
term can be done analytically and since Pb{^) is odd the 
Bose term does not contribute to the integral. |lj The 
net result is that 



5^= In 2 



dLupp{uj) 



(23) 



Additional information such as this may be imposed by 
modifying the prior or the likelihood function. Given the 
similarity of Eq. |2^ and Eq. ^0|, which appears as part 
of the likelihood function, we choose the latter approach. 
We introduce Soo as an additional datum with a relative 
error estimate as/Sao chosen to be approximately equal 
to the smallest relative error estimate of the energy data. 
We discuss our reasons for choosing this Ansatz in the 
appendix. 

With this Ansatz, we are in a position to employ the 
MEM. We write Kp from Eq. g as a linear combination 
of FpF + BpB- Eq. |l^ for x^ is modified similarly. We 
pick the default model in the following manner. We note 
first that it must be positive definite and integrable. We 
employ a Gaussian default model, which satisfies these 
criteria. 

Once the default model is selected, the method de- 
scribed in Ref. j|] may be applied straight away. To 
further illustrate the method, we now apply it to a non- 
trivial model. 



VI. EXAMPLE: 3D PERIODIC ANDERSON MODEL 

The 3d Periodic Anderson Model (PAM) is often used 
to investigate f-electron systems, where electronic corre- 
lations are important for the phenomena under study. It 
is a simplified lattice model in which the Coulomb in- 
teraction is limited in range to on-site interactions only 
and then only within one of the two bands. Neverthe- 
less, it is a rich model in which the interplay between 
delocalization (kinetic energy). Coulomb repulsion, Pauli 
exclusion, temperature, and electron density give rise to 
a wide variety of phenomena. 

The periodic Anderson Hamiltonian is 

ka k(T 

+UfY^{n^f^ - -){nify - -) 



+ '^efnif„ - fi'^{nifa + nida) ■ (24) 

We choose a simple cubic structure for which, 

efc = —2tdd [cos k^a + cos kya + cos k^a] , 

Vfc = —'2-tfd [cos /ca;a + cos kya + cos fczo] , (25) 

where a is the lattice constant. The dispersion of V}~ 
reflects our choice of near-neighbor (as opposed to on- 
site) hybridization of the / and d electrons. With on- 
site hybridization, the PAM is an insulator at half-filling, 
whereas with our intersite hybridization choice the half- 
filled, symmetric PAM is metallic. 




FIG. 1. Energy vs. temperature for various / — d hybridiza- 
tions tfd- The symbols mark QMC data, with errorbars. The 
solid lines mark the energy obtained via the MEM. There is 
an excellent agreement between the data and the MEM re- 
sults throughout the range of temperatures. This is further 
illustrated in the inset, where a magnified view of the lowest 
QMC temperature values for tfd = 0.2 is shown. 

The parameter values and the temperature T in this 
work are given in units of tdd- We take Uf = 6 and ex- 
plore a range of tfd and T values. QMC results for this 
model were obtained using the determinant algorithm, 
||12| which provides an exact treatment (to within sta- 
tistical errors and finite size effects) of the correlations. 
We further choose the symmetric PAM {fi = £f = 0, and 
thus half-filling: (riif) = (riid) = 1) in order to eliminate 
the QMC "sign problem," allowing accurate simulations 
at low temperatures. 

This version of the PAM has been studied in this pa- 
rameter regime and is known to undergo a sharp finite- 
temperature crossover at finite tfd ~ 0.6 — 0.8 with an 
associated, abrupt change in the free energy which is re- 
flected in the specific heat. ^^ These thermodynamic 
anomalies are believed to be signals of a zero-temperature 
metal-insulator phase transition in the f-band which is 



also seen as a finite-temperature crossover to a localized 
f-electron system. |l^,|lj] We will test our new technique 
by using it to reproduce the published work. 

Figure^ shows the energy obtained via the MEM us- 
ing Eq. |0 for various hybridizations tfd = 0.2,0.6,1.0 
along with QMC data and errorbars on the QMC data. 
There is an excellent agreement between the QMC data 
and the MEM results throughout the range of tempera- 
tures simulated by QMC. The quality of the agreement 
is further illustrated by the inset in Fig. n^, which shows 
a magnified view of the lowest QMC temperature values 
for tfd == 0.2. 




FIG. 2. Specific heat C{T), specific heat divided by tem- 
perature C{T)/T, and entropy S(T) (top figure), and energy 
E(T) (bottom figure) as a function of temperature T for a 
fixed hybridization tfd = 1.0. Representative error bars are 
shown. The dashed vertical line at T = 0.08 corresponds to 
the lowest QMC data point. A peak in C{T)/T appears at 
T ~ 0.15 due to singlet formation. At a higher tempera- 
ture, a hump appears in C{T)/T which is believed to be due 
to the suppression of charge fluctuations. The total entropy 
3.375 In 2 in the system obtained by integrating the image 
from the MEM (see the text) matches both the total entropy 
obtained from integrating C{T)/T (Eq. 0) and the value of 
Soo = 4 In 2 — ^o known for the model. 



Once one obtains the image, the specific heat C{T) is 
obtained by differentiating as in Eq. |^. Fig. g shows the 
specific heat C{T), the specific heat divided by the tem- 
perature C{T)/T, and the entropy S{T) in the top figure 
as a function of temperature T for a flxed hybridization 
tfd — 1.0. For comparison, the bottom figure shows the 
energy E(T). At this hybridization, singlets are known 
to form at low temperatures. [Q||] This is reflected in a 
peak in C{T)/T appearing at T « 0.15. The singlet for- 
mation peak is also visible in C(T). A smaller peak in 
C{T) at higher temperatures is believed to be due to the 
suppression of charge fluctuations in the /— band. 



The entropy may be found in two ways. First, S{T) 
may be obtained by integrating C{T)/T according to 
Eq. ||. This was done and is plotted in Fig. g. The 
entropy found in this manner saturates at high temper- 
atures at the infinite-temperature hmit Soo = 4 In 2 — 5'o 
known for this model. O] Second, the infinite tempera- 
ture entropy may calculated by integrating pp, Eq. |2^. 
The latter estimate is Soo = 3.375 In 2, which is also that 
obtained from integrating C{T)/T (Eq. |4|) and the value 
of 5*00 = 41n2 — 5o known for the model. 
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FIG. 3. Specific heat divided by temperature C{T)/T for 
various f ~ d hybridizations tfd, (a) tfd ~ 1.0, (b) tfd = 0.6, 
and (c) tfd = 0.2. The corresponding entropies from inte- 
grating C{T)/T (Eq. m are shown in panel (d). The dashed 
vertical line at T = 0.08 corresponds to the lowest QMC data 
point. At large hybridizations, tfd = 1.0, essentially all of the 
entropy in the system is quenched within the temperatures ac- 
cessed by the simulations (a). As the hybridization decreases, 
antiferromagnetic ordering at temperatures below those ac- 
cessed by the simulations becomes important and a substan- 
tial amount of entropy is not quenched within the simulated 
temperatures. The sum rule enforces the totaf entropy in the 
system, which appears in the resuits below the accessed tem- 
peratures, yet the most probable form of the specific heat for 
this regime cannot be precisely determined. This is reflected in 
large errorbars seen in panels (b) and (c) for the specific heat 
at temperatures below those accessed by the simulations. 

In addition to singlet formation at low temperatures 
for relatively large hybridizations tfd, the metallic PAM 
develops antiferromagnetic long-range-order (AFLRO) at 
low temperatures for small tfd- Hence, if one exam- 
ines C{T)/T for decreasing hybridization tfd, the singlet- 
formation peak should eventually disappear and a low- 
temperature peak corresponding to AFLRO should ap- 
pear in C{T)/T. Previous work has observed the disap- 
pearance of the singlet-formation peak, but did not access 
a sufficiently low temperature to observe the appearance 
of the AFLRO peak, pi] 



Figure || shows the specific heat divided by temper- 
ature C{T)/T for various f — d hybridizations tfd, (a) 
tfd ^ 1.0, (b) tfd = 0.6, and (c) tfd = 0.2. The corre- 
sponding entropies from integrating C{T)/T (Eq. Eh are 
shown in panel (d). Here, we observe the disappearance 
of the singlet peak with decreasing hybridization tfd- At 
tfd — 1.0 there is a substantial singlet formation peak. 
At tfd = 0.2 there is no singlet formation for any tem- 
perature accessed by the simulations, and possibly no 
singlet formation even at T = 0. p3| The intermediate 
hybridization tfd = 0.6 corresponds to a regime where 
singlet formation occurs suddenly for low temperatures, 
Q as is reflected in the shift of the singlet peak to lower 
temperatures in Fig g. 

The sum rule |l|l 



^oo=41n2-S'o 



(26) 



for the entropy enforces the total entropy in the system. 
This is an important feature of the method and satis- 
fying the sum rule is one check on whether the specific 
heat is physically reasonable. However, when the system 
does not quench all of the entropy by the lowest temper- 
ature accessed by the QMC simulation, one may worry 
whether enforcing the sum rule will push spurious en- 
tropy into the specific heat. This does not happen, as 
shown in Fig. H(d). Instead, this entropy beyond that 
quenched within the accessed temperatures goes below 
the lowest QMC temperature, where indeed it should go 
on physical grounds. However, then the extrapolation 
below the lowest QMC data point is unreliable. This un- 
reliability of the extrapolation to a regime where the sum 
rule has forced entropy below the QMC data is seen by 
the large error bars on the specific heat in this extrap- 
olation regime. That is, the method puts the entropy 
where it belongs, but then informs one that the results 
of the MEM in this regime are totally uncertain. This is 
an extremely desirable result. |16[ 



VII. SUMMARY 

We have described a novel technique to obtain the in- 
ternal energy as a function of temperature, as well as the 
specific heat, the entropy, and the free energy of a sys- 
tem using QMC energy data sampled at a small, finite 
set of temperature values. Our technique relies on proba- 
bility theory to obtain the most probable thermodynamic 
functions given the sampled QMC energy. The question 
of determining the number of fitting parameters, which 
plagues the other methods, is thereby removed. An en- 
tropy sum rule or other appropriate a priori information 
may also be used, if known. The technique was illustrated 
by applying it to the 3d Periodic Anderson Model. An 
important benefit of the technique is that it returns not 



only the thermodynamic functions, but also their uncer- 
tainties. This is a significant improvement over the prior 
techniques. 
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APPENDIX A: THE FORM OF THE ENERGY 

ANSATZ 

In section M we wrote the Ansatz for the energy as 

/oo 
dwwF(/3,w)/9F(w) 
-OO 

duJUjB{f3,uj)pB{cj), (Al) 



Fluctuation-Dissipation theorem. Therefore, we may re- 



where F and B are the Fermi- and Bose-distribution 
functions, respectively 



1 






1 



,/3w ■ 



(A2) 



Each integral in Al is a Fredholm integral equation of 



the first kind for the corresponding density with either 
F or B as the kernel. Since these kernels are continuous. 



the problem of inverting Eq. Al to find E{T) is ill-posed 
|l8| and must be regularized. MEM provides an efficient 
regularization method. In this appendix we motivate our 
Ansatz and discuss other possible functional forms. 

The Ansatz given in Eq. Al is not the only possible, or 
even sensible choice. In general one wants to choose as an 
Ansatz the one that captures the underlying physics of 
the problem. For models such as the PAM, we expect a 
density of Fermionic excitations represented by pp associ- 
ated with quasiparticle excitations. Since quasiparticles 
are conserved, their density should be positive definite 
which is required when the MEM formalism is employed. 

In addition, we expect a density pB of Bosonic exci- 
tations associated with collective behavior such as spin 
waves. With no Bosonic operators in the Hamiltonian, 
these Bosons are not conserved and have zero chemical 
potential. Thus, pb{uj > 0) > corresponding to the 
creation of such excitations, and pb{^ < 0) < corre- 
sponding to their destruction. This choice oi pp, pB, and 
K constrain the specific heat C{T) to be positive definite, 
C{T) > 0. Furthermore, PBi^^) is odd as required by the 



duce the second integral in Eq. Al to the range (0, oo) 
where pb(j-^) is positive semi-definite. 

An Ansatz should also provide a faithful representa- 
tion of the E{T) data. To explore this question, we note 
that for a finite-sized system E{T) is an analytic func- 
tion that can be expanded in a Taylor series in T around 
T = 0. Expanding F{P, uj) in a Sommerfeld low temper- 
ature expansion |l^ yields only even powers of T in the 
energy. In order to get odd powers of T that complete 
the Taylor series, the Bosonic kernel is required. Having 
the Taylor expansion for E(T), the remaining question is 
the positive-definite nature of the image p(w). Even the 
fact that the energy is monotonic does not yield a math- 
ematical constraint that pp and p_b(w > 0) are positive 
definite, yet in practice this constraint imposed by the 
MEM is not a limitation. 

It is possible to dispose of the Ansatz and use a general 
form for the energy. A faithful representation for systems 
in thermal equilibrium with a heat bath is provided by 



E{T) = - 



I^oo d^^p{^) exp{-uj/T) 
iZo dujp{uj) exp{-uj/T) 



(A3) 



where p{uj) > is the density of eigenstates of the Hamil- 
tonian. However, since the relationship between the den- 
sity and E is non-linear, the MEM algorithm described 
in Sec. I^ cannot be applied without significant mod- 
ification. Furthermore, the representation provided by 
Eq. |Al| in each case we have tested has provided a fit to 
within the measured error. Thus, the additional compli- 



cations associated with the use of Eq. A3 seem unneces- 
sary for this and similar systems. However, this general 
technique would allow extension of the method to Clas- 
sical Monte Carlo simulations and is something we are 
currently exploring. 
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